function [x0,place_dist,place_rad,dist0,prim_list,meta_list,inferred_data,size_projections] = getOriginSizes_noShed_v1_0(key,data,prim_acceptance_thresh)
%% getOriginSizes_noShed_v1_0
%  version 1.0
%  Author: Adeyinka Lesi
%  Date: 8/06/20
% get TIME_ZERO_SIZES using RK4 (step back in time using given velocities)
% need to handle tumors that are known to be present but are below
% detection limit. That implies that the specified size in x_init is
% irrelevant for these 'invisible' tumors and they will be placed based on
% a putative distribution
%% Version History
%  1.0: from getOriginSizes_v1_0,removing shedding from size trajectory
%  calculation

if(isfield(key,'QUIET_MODE') && ~key.QUIET_MODE)
    display(['> ' mfilename]);
end

% create key with zero shedding parameter
changes = struct('SHED_PARAMETER1',0,'SHED_PARAMETER2',0,...
    'CTC_PARAMETER1',0,'SHED_META_RATIO',1e3,'CTC_PARAMETER2',0,...
    'ORIGIN_TIME_CALCULATOR',@getOriginSizes_v1_0);
withShed_key = key;
key = changeKey_v3_1(withShed_key,changes);

x0 = key.TIME_ZERO_SIZES;
place_dist = key.PLACEMENT_DIST;
place_rad = key.PLACEMENT_RADIUS;
dist0 = key.TIME_ZERO_DIST;
prim_list = key.TIME_ZERO_ORIGIN;
meta_list = key.META_ORIGIN;
inferred_data = key.inferred_data;
size_projections = key.SIZE_PROJECTIONS;

    
% % initialize
% if(~exist('prim_acceptance_thresh','var'))
%     prim_acceptance_thresh = 3;
% end
% 
% nt = length(data.t);
% 
% if(isfield(data,'cdf_adj'))
%     p_adj = data.cdf_adj(1);
% else
%     p_adj = 0;
% end
% x_t = cell(1,nt);
% dist = cell(1,nt);
% cdf = cell(1,nt);
% x_t{1} = data.x(data.selector{1});
% dist{1} = data.dist(data.selector{1},1);
% cdf{1} = p_adj+data.cdf(end,1)-[0; data.cdf(data.selector{1}(1:end-1),1)];
% x0_size = length(x_t{1});
% for n = 2:nt
%     if(isfield(data,'cdf_adj'))
%         p_adj = data.cdf_adj(n);
%     else
%         p_adj = 0;
%     end
%     x_t{n} = data.x(data.selector{n});
%     dist{n} = data.dist(data.selector{n},n);
%     cdf{n} = p_adj+data.cdf(end,n)-[0; data.cdf(data.selector{n}(1:end-1),n)];
%     x0_size = max(x0_size,length(x_t{n}));
% end
% 
% inferred_data = data;
% 
% temp_data = struct('x',cell(1,nt),'dist',cell(1,nt),'ccdf',...
%     cell(1,nt),'step_ccdf',cell(1,nt),'labels',cell(1,nt),...
%     'detect_time',cell(1,nt));
% 
% x0 = zeros(x0_size,1);
% dist0 = zeros(x0_size,1);
% prim_list = zeros(x0_size,3);
% known_list = zeros(x0_size,5);
% meta_list = zeros(x0_size,5);
% backtrack_sizes = cell(1,nt);
% origin_times = cell(1,nt);
% nx0 = 0;
% nkwn = 0;
% nmet = 0;
% place_dist = zeros(size(x0));
% place_rad = zeros(size(x0));
% i_selected = cell(1,nt);
% x_t1_low = 0;
% ntum_inc = 0;
% low_check1 = 0; % low limit for first time point
% low_check2 = 1; % low limit for other time points
% for i = 1:nt
%     x_ti = x_t{i};
%     x0i = zeros(length(x_ti),1);
%     place_disti = zeros(size(x0));
%     place_radi = zeros(size(x0));
%     ts = data.t(i);
%     prev_sizes_store = zeros(length(x_ti),i);
%     birth_times = zeros(length(x_ti),1);
%     if(ts>0)
%         dt = ts/round(ts/key.TIME_STEP)/2;
%         for j = 1:length(x_ti)
%             if(x_ti(j)==0)
%                 % indication disappeared tumor
%                 x0i(j) = 0;
%                 prev_sizes = zeros(1,i);
%                 prev_sizes_store(j,:) = prev_sizes;
%                 birth_times(j) = -1;
%             else
%                 [traj,ttraj] = getSizeTrajectory_AlwaysSmaller_v4_0(key,x_ti(j),-dt,ts);
%                 if(traj(end) < x_ti(j))
%                     x0i(j) = max(1,traj(end));
%                 else
%                     % net negative velocity
%                     % 1) parameters don't make sense
%                     % 2) will set initial point to 1 (for now) - effect of this
%                     % depends on low_checks
%                     x0i(j) = 1;
%                 end
%                 % store predicted sizes at previous times
%                 prev_ti = length(traj)-round(data.t(1:i)/dt);
%                 prev_sizes = traj(prev_ti);
%                 prev_sizes(prev_sizes<1) = 0;
%                 prev_sizes_store(j,:) = prev_sizes;
%                 % get birth time
%                 birth_ind = find(traj>1,1,'last');
%                 if(~isempty(birth_ind))
%                     birth_times(j) = ttraj(birth_ind);
%                 end
%             end
%         end
%     else
%         for j = 1:length(x_ti)
%             x0i(j) = x_ti(j);
%             prev_sizes = x_ti(j);
%             prev_sizes_store(j,:) = prev_sizes;
%             %             if(x_ti(j) >= lowx)
%             %                 x0i(j) = x_ti(j);
%             %             else
%             %                 x0i(j) = 1+0.5*(lowx-1);
%             %                 place_disti(j) = unk_dist_key.(unk_dist);
%             %                 place_radi(j) = 0.5*(lowx-1);
%             %             end
%         end
%     end
%     
%     % test to see what initial tumors should be accepted
%     
%     % get rid of x0 values < low_check1
%     i_above = find(x0i>low_check1);
%     if(isempty(i_above))
%         fprintf('Using 0 initial tumor values from time %i\n',i);
%     else
%         if(i == 1)
%             if(length(i_above) < length(x0i))
%                 fprintf('Using %i out of %i initial tumor values from time %i\n',length(i_above),length(x0i),i);
%             else
%                 fprintf('Using all %i initial tumor values from time %i\n',length(x0i),i);
%             end
%             
%             nai = length(i_above);
%             x0(nx0+1:nx0+nai) = x0i(i_above);
%             dist0(nx0+1:nx0+nai) = dist{i}(i_above);
%             prim_list(nx0+1:nx0+nai,:) = [ones(nai,1)*i i_above x_ti(i_above)];
%             place_dist(nx0+1:nx0+nai) = place_disti(i_above);
%             place_rad(nx0+1:nx0+nai) = place_radi(i_above);
%             i_selected{i} = i_above;
%             nx0 = nx0+nai;
%             
%             % this is not needed...
%             %             known_sizes = x0i(i_above);
%             %             known_origin_time = birth_times(i_above);
%             %             known_list(nkwn+1:nkwn+nai,:) = [ones(nai,1)*i i_above known_sizes known_origin_time];
%             %             nkwn = nkwn+i_above;
%             
%             if(nx0 > 0)
%                 x_t1_low = x0(nx0); % lowest size at time 1
%                 ntum_inc = sum(dist{1}(i_above)); % normalized number of tumors included
%             elseif(~isempty(x_t{1}))
%                 x_t1_low = x_t{1}(end);
%                 ntum_inc = 0;
%             else
%                 x_t1_low = inf;
%                 ntum_inc = 0;
%             end
%             %             x_t1_low = inf;
%             %             n_t1 = 0;
%         else
%             % want to select init values larger than 1 for the case when
%             % there are excess tumors at later times
%             
%             if(i <= prim_acceptance_thresh)
%                 prim_sel = find((x0i>low_check2).*(cdf{i}>ntum_inc)); % labeled as new (formerly undetected) primary
%                 meta_sel = find((x0i>0).*(x0i<=low_check2).*(cdf{i}>ntum_inc)); % this data is attributed to new (formerly undetected) metastases
%                 known_sel = find((x0i>0).*(cdf{i}<=ntum_inc)); % this data is attributed to previously identified tumors (prim or met)
%             else
%                 prim_sel = [];
%                 meta_sel = find((x0i>0).*(cdf{i}>ntum_inc)); % this data is attributed to new (formerly undetected) metastases
%                 known_sel = find((x0i>0).*(cdf{i}<=ntum_inc)); % this data is attributed to previously identified tumors (prim or met)
%             end
%             if(length(prim_sel) < length(x0i))
%                 fprintf('Using %i out of %i initial tumor values from time %i\n',length(prim_sel),length(x0i),i);
%             else
%                 fprintf('Using all %i initial tumor values from time %i\n',length(x0i),i);
%             end
%             
%             % update primaries
%             nai = length(prim_sel);
%             x0(nx0+1:nx0+nai) = x0i(prim_sel);
%             dist0(nx0+1:nx0+nai) = dist{i}(prim_sel);
%             prim_sizes = x_ti(prim_sel);
%             prim_list(nx0+1:nx0+nai,:) = [ones(nai,1)*i prim_sel prim_sizes];
%             place_dist(nx0+1:nx0+nai) = place_disti(prim_sel);
%             place_rad(nx0+1:nx0+nai) = place_radi(prim_sel);
%             i_selected{i} = prim_sel;
%             nx0 = nx0+nai;
%             
%             % update known primaries and metastasis
%             nki = length(known_sel);
%             known_sizes = x_ti(known_sel);
%             known_dist = dist{i}(known_sel);
%             known_origin_time = birth_times(known_sel);
%             known_list(nkwn+1:nkwn+nki,:) = [ones(nki,1)*i known_sel known_sizes known_dist known_origin_time];
%             nkwn = nkwn+nki;
%             
%             % update metastases
%             nmi = length(meta_sel);
%             meta_sizes = x_ti(meta_sel);
%             meta_dist = dist{i}(meta_sel);
%             meta_origin_time = birth_times(meta_sel);
%             meta_list(nmet+1:nmet+nmi,:) = [ones(nmi,1)*i meta_sel meta_sizes meta_dist meta_origin_time];
%             nmet = nmet+nmi;
%             
%             % update total number of tumors in distribution
%             ntum_inc = ntum_inc+sum(dist{i}(prim_sel))+sum(dist{i}(meta_sel));
%         end
%     end
%     backtrack_sizes{i} = prev_sizes_store;
%     origin_times{i} = birth_times;
% end
% 
% if(nx0 == 0)
%     x0(1) = 1;
%     nx0 = nx0+1;
%     warning('None of the initial tumor values were used. Using one size 1 tumor');
% end
% 
% x0 = x0(1:nx0);
% [~,ord] = sort(x0(1:nx0));
% x0 = x0(ord);
% dist0 = dist0(ord);
% prim_list = prim_list(ord,:);
% place_dist = place_dist(ord);
% place_rad = place_rad(ord);
% known_list = known_list(1:nkwn,:);
% meta_list = meta_list(1:nmet,:);
% 
% size_projections = struct('backtrack_sizes',backtrack_sizes,'origin_times',origin_times);
% 
% % construct new data distribution (with inferred primaries and with and 
% % without inferred meta
% % add backtrack_sizes and origin_times to structure
% 
% % construct distributions
% % 1) known (visible prim and meta)+inferred prim
% % 2) known (visible prim and meta)+inferred prim+inferred meta
% nxt = 0;
% for i = 1:nt
%     known_i = known_list(:,1)==i;
%     prim_i = prim_list(:,1)==i;
%     meta_i = meta_list(:,1)==i;
%     infp_i = prim_list(:,1)>i;
%     infm_i = (meta_list(:,1)>i)&(meta_list(:,end)<=data.t(i));
%     n_known = sum(known_i);
%     n_newprim = sum(prim_i);
%     n_newmeta = sum(meta_i);
%     n_infprim = sum(infp_i);
%     n_infmeta = sum(infm_i);
%     
%     xi_lab = [ones(n_infmeta,1)*2; ones(n_infprim,1)*1; zeros(n_newmeta+n_newprim+n_known,1)];
%     xi_orig = zeros(size(xi_lab));
%     xi = zeros(size(xi_lab));
%     disti = zeros(size(xi_lab));
%     nxi = 0;
%     xi(end-nxi-n_known+1:end-nxi) = backtrack_sizes{i}(known_list(known_i,2),i);
%     disti(end-nxi-n_known+1:end-nxi) = dist{i}(known_list(known_i,2),1);
%     nxi = nxi+n_known;
%     xi(end-nxi-n_newprim+1:end-nxi) = backtrack_sizes{i}(prim_list(prim_i,2),i);
%     disti(end-nxi-n_newprim+1:end-nxi) = dist{i}(prim_list(prim_i,2),1);
%     nxi = nxi+n_newprim;
%     xi(end-nxi-n_newmeta+1:end-nxi) = backtrack_sizes{i}(meta_list(meta_i,2),i);
%     disti(end-nxi-n_newmeta+1:end-nxi) = dist{i}(meta_list(meta_i,2),1);
%     nxi = nxi+n_newmeta;
%     xi_orig(end-nxi+1:end) = ones(nxi,1)*i;
%     
%     infp_inds = find(infp_i);
%     for j = 1:n_infprim
%         xi(end-nxi) = backtrack_sizes{prim_list(infp_inds(j),1)}(prim_list(infp_inds(j),2),i);
%         disti(end-nxi) = dist{prim_list(infp_inds(j),1)}(prim_list(infp_inds(j),2),1);
%         xi_orig(end-nxi) = prim_list(infp_inds(j),1);
%         nxi = nxi+1;
%     end
%     infm_inds = find(infm_i);
%     for j = 1:n_infmeta
%         xi(end-nxi) = backtrack_sizes{meta_list(infm_inds(j),1)}(meta_list(infm_inds(j),2),i);
%         disti(end-nxi) = dist{meta_list(infm_inds(j),1)}(meta_list(infm_inds(j),2),1);
%         xi_orig(end-nxi) = meta_list(infm_inds(j),1);
%         nxi = nxi+1;
%     end
%     
%     % sort by size
%     [~,ord] = sort(xi);
%     xi = xi(ord);
%     disti = disti(ord);
%     xi_lab = xi_lab(ord);
%     xi_orig = xi_orig(ord);
%     
%     % calculate cdf...
%     if(isfield(data,'cdf_adj'))
%         p_adj = data.cdf_adj(i);
%     else
%         p_adj = 0;
%     end
%     cdfi = getStepIntegral(xi,disti);
%     ccdfi = p_adj+cdfi(end)-[0; cdfi];
%     ccdf2i = zeros(2*length(xi),3); % double cdf for making plots
%     ccdf2i(2*(1:length(xi))-1,1) = xi;
%     ccdf2i(2*(1:length(xi)),1) = xi;
%     ccdf2i(2*(1:length(xi))-1,2) = ccdfi(1:end-1);
%     ccdf2i(2*(1:length(xi)),2) = ccdfi(2:end);
%     ccdf2i(2*(1:length(xi))-1,3) = xi_lab;
%     ccdf2i(2*(1:length(xi)),3) = xi_lab;
%     
%     temp_data(i).x = xi;
%     temp_data(i).dist = disti;
%     temp_data(i).ccdf = ccdfi(1:end-1);
%     temp_data(i).step_ccdf = ccdf2i;
%     temp_data(i).labels = xi_lab;
%     temp_data(i).detect_time = xi_orig;
%     
%     nxt = nxt+nxi;
%     
%     %     % plott
%     %     figure;
%     %     semilogx(x_t{i},cdf{i},'--ob',ccdf2i(:,1),ccdf2i(:,2),'-g',...
%     %         ccdf2i(ccdf2i(:,3)==1,1),ccdf2i(ccdf2i(:,3)==1,2),'sr',...
%     %         ccdf2i(ccdf2i(:,3)==2,1),ccdf2i(ccdf2i(:,3)==2,2),'dr');
%     %     if(n_infprim>0 && n_infmeta>0)
%     %         legend('Original','New','Inferred Primaries','Inferred Metastasis');
%     %     elseif(n_infprim>0)
%     %         legend('Original','New','Inferred Primaries');
%     %     elseif(n_infmeta>0)
%     %         legend('Original','New','Inferred Metastasis');
%     %     else
%     %         legend('Original','New');
%     %     end
%     %     title(sprintf('Day %i',data.t(i)));
%     %     axis([1 3e7 0 2]);
% end
% 
% % convert inferred distributions to same format as 'data' structure
% % dist: [263�7 double]
% % a: [263�1 double]
% % t: [1 3 5 7 9 11 13]
% % x: [263�1 double]
% % cdf: [263�7 double]
% % selector: {[48�1 double]  [57�1 double]  [50�1 double]  [35�1 double]  [24�1 double]  [19�1 double]  [18�1 double]}
% % max: 2.8861e+07
% % lowa: 0.4393
% % lowx: 1.1713
% % conv: [1�1 struct]
% % xlabel: 'Tumor Size'
% % ylabel: 'Number of Tumors per Body'
% % plot_indices: [1 2 3 4 5 6 7]
% % weights: [1 1 1 1 1 1 1]
% xd = zeros(nxt,1);
% % ad = zeros(nxt,1);
% distd = zeros(nxt,nt);
% cdfd = zeros(nxt,nt);
% labelsd = zeros(nxt,nt);
% selectord = cell(1,nt);
% ninc = 0;
% for i = 1:nt
%     xd(ninc+1:ninc+length(temp_data(i).x)) = temp_data(i).x;
%     distd(ninc+1:ninc+length(temp_data(i).x),i) = temp_data(i).dist;
%     labelsd(ninc+1:ninc+length(temp_data(i).x),i) = temp_data(i).labels;
%     ninc = ninc+length(temp_data(i).x);
% end
% % sort
% [~,ord2] = sort(xd);
% xd = xd(ord2);
% ad = data.conv.v2a(xd);
% distd = distd(ord2,:);
% labelsd = labelsd(ord2,:);
% % get selector and cdf
% for i = 1:nt
%     cdfd(:,i) = getStepIntegral(ad,distd(:,i));
%     selectord{i} = find(distd(:,i)>0);
% end
% % place correct values into structure
% inferred_data.dist = distd;
% inferred_data.a = ad;
% inferred_data.x = xd;
% inferred_data.cdf = cdfd;
% inferred_data.selector = selectord;
% inferred_data.max = max(xd);
% % store new format data into output
% inferred_data.labels = labelsd;
% inferred_data.alternate = temp_data;

end
